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We present an accurate and rapid solution of Poisson's equation for space-filling, arbitrarily- 
shaped, convex Voronoi polyhedra (VP); the method is O(Nvp), where Nvp is the number of 
distinct VP representing the system. In effect, we resolve the longstanding problem of fast but 
accurate numerical solution of the near-field corrections (NFC), contributions to each VP potential 
from nearby VP - typically involving multipole-type conditionally-convergent sums, or fast Fourier 
transforms. Our method avoids all ill-convergent sums, is simple, accurate, efficient, and works gen- 
erally, i.e., for periodic solids, molecules, or systems with disorder or imperfections. We demonstrate 
the method's practicality by numerical calculations compared to exactly solvable models. 

PACS numbers: 41.20.Cv, 71.15.Dx 



I. INTRODUCTION 



Poisson's equation describes the electrostatics by re- 
lating a charge distribution to the potential contingent 
upon the boundary conditions. An accurate solution of 
Poisson's equation is critical in various areas of chemistry 
and condensed-matter physics. In ah initio electronic- 
structure methods, the Poisson equation is solved re- 
peatedly, and concomitantly parallel to the Schrodinger's 
equation. As such, computational time for solving Pois- 
son equation is always a concern. Although a number of 
proposals exist, most suffer from shortcomings that affect 
accuracy and speed, and the ability to scale to large sys- 
tem sizes efficiently. Here we provide an exact treatment 
of Poisson's equation and its accurate and efficient nu- 
merical solution of the potential and Coulomb energy of 
systems described by arbitrarily-shaped, convex, space- 
filling VP in any site-centered method. Our new ap- 
proach scales linearly with the number of VP, and avoids 
mathematical and numerical issues associated with pre- 
vious methods, particularly multipole approaches. In 
historical context, we provide an efficient and accurate 
means to compute the so-called "near-field corrections" 
(NFC), a problem not fully resolved so far. 

Typically, the electrostatic potential at a point in a 
convex VP is given by two contributions,^ ^ namely, (i) 
an intracell potential arising from the charge density 
within a VP (p*^°^ in rig) and (ii) an intercell potential 
arising from all other p^^' in il^'s, see Fig. [I] In general. 
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where p^^^ is a truncated density centered at site R. 
Computational time in most methodJ^^ arise from the 
use of L = {/,to} multipole (spherical-harmonics Yl(y)) 
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FIG. 1. (Color online) Two VPs (fio and S7_r) separated by 
vector R with overlapping bounding spheres with radius fbs- 
For rmin < r < rss, NFC are needed, tmt is the inscribed 
sphere radius (not drawn for clarity). 



expansions. Evaluation of intercell potential (term two 
in Eq. ([T])) is the most tricky, and our main focus. Of- 
ten, as a first step, the Green's function |r — (r' + R)|^^ 
is expanded in l^'s in terms of r< (e.g., |r|) and 
Ir'-t-Rr 



(e.g. 



see Sec. 
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attempting to sep- 
arate two of three {r,r',R) degrees of freedom. In 
most existing methods,'^lt2lan additional multipole expan- 
sion of (r' + R) is performed yielding conditionally- 
convergent nested L-sums (internal vs. external: /JJJ*^ > 
3^5jf*x) due to the nearest-neighbor sites, and relevant in 
the light shaded (pink) region in Fig. [T] Such nested 
sums are numerically expensive and ill convergent, even 
more so for distorted (asymmetric) cells. Numerical inef- 
ficiency also arises from any use of VP shape functions]^ 
which utilize Y^'s to expand VP shapes to facilitate VP 
integrations; again, these are costly (and inaccurate) due 
to the large L-sums (^'"'>>3ZJ^ax) required. For "muffin- 
tin" potentials varying only inside ta/t (Fig. [I]), these 
issues are moot as no conditional expansions are needed; 
the "atomic sphere approximation" ignores these errors. 



2 



Thus, for arbitrarily-shaped, convex, space-fiUing VP, 
we derive the set of integral equations that permit us 
to eliminate all previous computational bottlenecks and 
convergence issues to solve Poisson's equation by employ- 
ing isoparametric integration,^'' valid for rapidly varying 
and/or decaying integrands, while providing a dramatic 
savings of computational time, e.g., 10^ in time and 
10^ in accuracy over the shape-functions! The method 
permits site-specific quantities to be calculated rapidly, 
scales linearly with the number of VP Nyp and is easily 
parallelized. Unlike the Full-potential Linear Augmented 
Plane- Wave (FLAPW) method. Fast Fourier Transforms 
(FFTs), which limit scaling to large systems, are not 
needed. To prove these points explicitly, we compute ex- 
ample integrals for potential a nd C oulomb energy from 
analytic charge- density models.'^n] 



cause it is faster albeit approximate.!^ 



III. A COMPUTATIONALLY EFFICIENT AND 
ACCURATE POISSON SOLVER 

A proposal by Nicholson and SheltorP is conceptually 
easy, although it suffers also from convergence issues - 
both multipoles and shape-functions. We use a key idea 
from their work but, uniquely in our derivation, avoid 
any expansions used in prior approaches, made possible 
by isoparametric integration.121 

To start, using L = {l,m} as a composite index, we 
express the solution of Poisson's equation a^ 

/max 

- E [^^(0 + «L r'] Y^i?), r < tbs (2) 



II. BACKGROUND 

To solve Poisson's equation for site-centered meth- 
ods, various techniques have been developed. Gonis et 
al^ introduced a technique (modified later by Vitos et 
al^ based on shifting (and back-shifting) the neighbor- 
ing cells by a vector b that eliminates the conditionally- 
convergent expansion related to these neighbors, but re- 
quires additional L sums; the technique converges very 
slowly versus Lmax because internal sums are large, e.g., 
^max > 3ZJ5fax; additionally, b is a parameter that must be 
chosen wisely and depends on crystal symmetry. Others^ 
used shape-functions making the VP integrations very 
fast for a Yi,-basis but the expansion is slowly conver- 
gent (i.e., > 30), with hmited accuracyf^Sl SchadleP 
proposed corrections to the usual multipole expansion via 
a conditionally-convergent formula due to Sackj^^ how- 
ever, these corrections do not satisfy Laplace's equation. 
Zhang et aLl^ converted VP integrals to surface integrals, 
avoiding most conditionally-convergent sums; however, 
it is not automated for complex geometries, and con- 
cerns remain about degeneracies for their set of linear 
equations. For FLAPW, WeineriP avoided these issues 
via Yt-basis in MT-spheres and interstitial plane-waves; 
however, to obtain a smooth density (for a chosen set of 
MT radii) a large number of plane waves (Npw>30,000) 
and Yj^'s {lmax > 8) are required, and one never obtains 
VP-specific properties. FFTs are then needed, scaling as 
2Npw^O(7(Npw), with specialized programming for large 
system sizes. For Linea r Combination of Atomic Orbital 
(LCAO) methods,'A^E3I various atomic bases (e.g., Gaus- 
sian orbitals) are used in different regions of space to 
study molecules and clusters. Gaussian-orbital methods 
do not necessarily require partitioning of space because 
Poisson's equation can be solved analytically (or in terms 
of incomplete Gamma functions) on any mesh of points. 
However, a significant advantage could be achieved by a 
method that solves Poisson's equation numerically and 
accurately; for example, some Gaussian-orbital codes re- 
sort to least-square fits to solve Poisson's equation be- 
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p'f^{r) is the extended charge density inside the circum- 
scribing (or bounding) sphere of radius tbs of the central 
cell Oq in Fig. [l] The radial function WL(r) is the con- 
tribution to the potential within a distance r from origin 
of f2oj which is given by 



dx 



„;+2 fix 
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and which is bounded, i.e., Wi(r — ;> 0) = 0, and finite for 
any r < rBs, and, therefore, easily integrated. 

The intracell potential is the first term in Eq. ^ , while 
the intercell potential was expressed as aLr'li to make 
apparent a mathematical "trick" (assignment of equality) 
used below. Here is an unknown coefficient depend- 
ing on the charge distribution of the system. The main 
objective is to determine ol, which, if known, would give 
the potential at any point inside the central rss sphere. 

The problem in calculating (r) directly in Eq. ^ 

is the need to assume (particularly for multipole ap- 
proaches) the geometric condition 

/ < K 



r < r 



■RL 
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which is not fulfill ed in the so-called moon region between 
the near VP cells ,'212121 shown by light (pink) shading in 
Fig.[lj or, in other words, the complement of the VP and 
its bounding sphere with radius rss- A cell centered at 
R is a near-cell of the central one if i? < + r'^gg ■ 
Incorrect contributions to the potential arise from near 
VP beyond a radius rmin, which have been often ignored 
or badly approximated. If, however, we limit ourselves 
to r < Tmin (Fig. [ij, the geometric condition Eq. Q 
is valid and the potential ([T]) can be calculated easily. 
The unknown coefficients aL can be then determined by 
equating Eqs. ^ and ^ within r^m- 

Now, following this line of reasoning, with r < r,„in 
r < |r' + R|, term two of Eq. ([T]) can be expressed a^ 

^aLr'lK?), where (5) 
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Rapidly varying and/or decaying integrand, as in Eq. 
Q , over general VP can be calculated accurately and fast 
with an isoparametric numerical quadrature methodl^l 
with a naly tically-known points and weight. (Other 
methods^ for performing integrals also works well, al- 
beit not as efficiently). A critical side point: no expan- 
sion (or FFT) of the integrand in Eq. (6]) is necessary, 
eliminating all previous computational bottlenecks and 
convergence issues. A rigorous example is provided in 
SecHVl 

Then, with p — )■ p'^^ for r < r-axia (the spherically sym- 
metric regime) , the first term of Eq. ([T]) is simplified as 
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Substituting Eqs. ([t]) and (|3| into Eq. ([T]) and comparing 
it with Eq. ^ yields for all rmin < r' < tbs (the 



remaining space), i.e., 
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Equation ([8]) is our central result. It serves to calculate 
accurately V^"*'^''(r) with the necessary NFC, given by the 
integral term. This NFC is non-zero only beyond r^^^ 
[p p'^^ for r < rmin) and pronounced in the "moon 
region" (r ^ flo and |r| < rss of fio)- 

Notably, knowing V^'^'^Xr < r^in) gives aL and, thus, 
V(r) everywhere in space via Eq. ([2]), which is ultimately 
the "trick" . Finally, the cell integrations in Eq. Q , which 
can exhibit rapidly varying and/or decaying integrands, 
needs to be performed by an accurate and fast integration 
method over arbitrarily-shaped VP, which is satisfied by 
a recently proposed isoparametric integration .ESI 

NFC provide the correct l/^"*'" (r) from the near-cells, 
and are the motivation behind previous methods.'^'SEHal 
Unlike existing schemes that address NFC, our deriva- 
tion is simple and provides an efficient, fast and accurate 
solution of Poisson's equation. 

In historical context, the ill-convergent sums in 
other methods arise from traditionally expanding 



le y£(r' + R)/|r' + R|'+i in Eq. i.e., for all r' < R, 
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which separate r' and R creating a multipole-type ex- 
pression via Eq. ^ with large internal, conditionally- 
convergent sums (L'"*). The convergence of such expan- 
sions (involving Gaunt coefficients C^^,,) is sensitive to 
the location of r' when R is a near-cell vector, being es- 
pecially difhcult to converge if r' lies, e.g., near one of 
the corners of the VP. To achieve a minimal level of con- 
vergence (e.g., 10"'*), the number of Us required is huge 
{I > 70) even for highly symmetric VP, such as fee and 
bcc! These errors are often ignored. 

For completeness, we note that the expansion neces- 
sary for the electrostatic potential for general charge dis- 
tributions in terms of spherical harmonics, like Eq. ^ 
has a long history which continues. For example, for one- 
and two-center Coulomb potentials, Buehler addressed 
spherical distributions,^ and Fontana addressed discrete 
distributions,-^ Jansen provided a tensor formalism for 
multipole expansionsj^^ however. Sack's results are well- 
known, as discussed in the Background section,^^ and 
often revisited*^ because of the use of hypergeomet- 
ric functions, which even Sack did later.ISS Nonetheless, 
all the results have extensive sums that are conditionally 
convergent. 

Finally, Gonis et aZ.'^Ei acknowledged that, in their 
method for solving Poisson equation, the /-convergence 



depends sensitively on the choice of the shifting vector b 
that mathematically moves the central site far enough 
away from the remaining nearest-neighbor sites such that 
the usual r< and r> spherical harmonic expansions are 
valid for all r within Oq; however, such a shifted expan- 
sion requires a very large intern al L s um for full conver- 
gence. In the resulting equation^^'^'^*' the shifting vector 
adds another conditionally-convergent summation, with 
multiply nested L sums. For large Ps, convergence fur- 
ther suffers due to the non-vanishing high multipole 
moments constructed from the shape function, giving 
slowly convergent inner sums for near cells and high Zout ■ 
Our method is free from such issues. 



IV. RESULTS AND DISCUSSION 

To illustrate the accuracy of our method, we present 
results for two distinctly different cases. First, an elec- 
tronic charge density model by van W. Morgan,^ in which 
all results can be derived and evaluated analytically, and 
which mirrors the collective densities of real atoms. Sec- 
ond, we address the well-known "Madelung" problem (a 
jellium-like model), which has a closed- form solution us- 
ing Ewald's method, but requires numerical evaluation 
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due to appearance of non-elementary special functions 
(error functions), as detailed over decades and presented 
in Slater's bool^^lfrom the work of Slater and de Cicco.'^S 



A. van Morgan density model 

To illustrate the accuracy of our method for the po- 
tential and Coulomb energy, we chose an analytic model 
by van W. Morgan,^ whose charge density is given by 



iT,,.r 



(10) 



B is an arbitrary constant (set to 1) and T„ (with mag- 
nitude |T|) are reciprocal-lattice vectors of the system 
under consideration, see Ref. for more details with 
the derived expression given in its appendix. The exact 
potential for such a charge distribution is 



V{r) = 47rp(r)|r| 



Vo, 



(11) 



where Vq is an arbitrary constant. Also, the Coulomb 
energy for VP unit-cell volume fio is 



U = 



1 



p{r)V{r)dr 



exact 27rr2o ■r-^ 



(12) 



n 



This charge-density model, which mimics real (collec- 
tive atomic-centered density) behavior provides a rigor- 
ous (exact) test, not possible in applications to a "real" 
system. 



For the density given by Eq. ( 10 1, we evaluate the first 
key integral quantity, provided in Eq. Q. Table |l] shows 
the coefficients (Eq. (|8|) with respect to the num- 
ber of Gauss points {Nq} to achieve 6 decimal place 



TABLE I. aim calculated via Eq. ([8| for fee {R is summed to 
8"* neighbor shell). {Ng} is the number of Gauss points per 
X, y, z direction for 6 decimal place accuracy, aoo does not 
match the exact result due to an overall constant of integra- 
tion, which depends on the crystal symmetry under consider- 
ation; however, it does not affect r-dependence. 
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FIG. 2. (Color online) V(r), relative to a constant, for various 
/max along high-symmetry directions in WS-cells of fee (top) 
and bcc (bottom) for van W. Morgan model. 



accuracy for various L = {Z,m}. The numerically cal- 
culated q;l are compared with the analytical exact ex- 
pression (right most column in Table |l]) given by, with 
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dm + VaSlo (13) 



and ji are the spherical Bessel function. In spite of the os- 
cillatory angular dependence in Eq. with Z-dependent 
spatial decay, the increase in Nq required with larger Vs is 
not significant, and, hence, the isoparametric integration 
method used remains fast. Only the aoo coefficient is not 
produced correctly, see Table|lj however, we note that (1) 
aoo is highly sensitive to the boundary conditions in the 
r — >■ cx) limit and how this limit is taken, see discussion 
by van W. Morgan ( appendix) or by Leeuw,'22! which 
nonetheless can be solved by standard Ewald techniques; 
and (2) the potential is defined up to an arbitrary con- 
stant generally, as used in most electronic-structure codes 
to advantage. Hence, the error in aoo does not impact 
the key spatial-dependence of the potential required. 

In Fig. [2] we compare V{y) calculated from Eq. ^ 
for ^max = 0,4,6,8, 10 with that of the exact result for 
fee and bcc lattices. The potential converges rapidly in /, 
with I — 8 results agreeing well with V^xact- The quality 
of agreement between the curves depends on the direction 
inside the VP cell, with ^-convergence slower for points 
near cell boundaries. For instance, H (P) symmetry point 
is the near (far) part of the fee VP, and X (L) is near (far) 
part of the bcc VP. Figure [3] shows the convergence of 
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the potential at these symmetry points versus Zmax; the 
potential at Z,„ax — 6 already converges within 0.1% of 
the exact result. Unlike previous approaches, our method 
requires just one converged L-sum (Zmax — 6 — 8), giving 
a significant speed up. 

The slower rate of Z-convergence near the cell boundary 
mainly arise due to larger NFC (integral term in Eq. ([s])) 
in this region, see Fig. [4j where the NFC to the potential 
for an fee lattice are shown along the two symmetry direc- 
tions with /max = 8. The potential within r^^^ with(out) 
NFC are the same as the exact result, as expected, and 
only beyond rmin does the correction grow. The NFC, 
although apparently small, are very important in get- 
ting the correct result, and are larger in less-symmetric 
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FIG. 4. (Color online) For fee, the potential with (without) 
NFC along F— H and F— P for van W. Morgan model. Solid 
eurves match with the exact results, vmt is an inscribed MT- 
sphere radius. 



structures, which may require a higher L-sum to con- 
verge. Moreover, the NFC for high L's are actually very 
large but compensated by the ql coefficients, and, at 
small L's the NFC are similar in magnitude to the ql's 
in most cases, making the integral term in Eq. (|8| critical 
to achieve the correct result. 

Figure [5] shows the convergence of Coulomb energy U 
versus Zmax for fee and bee lattices, compared to the van 
W. Morgan exact result. Without the NFC, the error is 
~10 mRy for fee and ~6 mRy for bee cases, and do not 
improve with higher L's. (No systematic error cancella- 
tion is possible, e.g., for Ufcc-Ubcc-) Unlike the potential, 
the Coulomb energy is almost exact by Imax — 6, because 
V — Vexact oscillates about zero for a given r as a func- 
tion of {9, (j)) and these contributions mostly cancel when 
integrated over the VP, which may be true for most cases. 



B. Madelung's Problem 

The Madelung "jellium" model consists of a constant 
electronic (negative) charge density throughout space, 
~Pa {Po = Z/Qq) which integrates to —Z, compensated 
by an ordered array of positive nuclear point charges +Z 
at atom-center positions R„, providing charge neutral- 
ity on average, locally (within a Voronoi or Wigner-Seitz 
cell) and globally. The total density then is 



ptot = ^ V'(5(r- R„) - Po. 



(14) 
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FIG. 6. (Color online) rV(r) for various /max along high- 
symmetry directions in WS-cells of fee (top) and bee (bottom) 
for Madelung jellium model. 



Via the Ewald methotpSi a compensating set of positive 
and negative Gaussian charge distributions are used, i.e., 



P?{r) 



7r3/2' 



(15) 



Similar to the van W. Morgan case, the accuracy of 
the potential for this jellium model varies along the high- 
symmetry directions, being worse at the H, P point for 
fee, and X, L point for bcc case, hence, requiring a higher 
L-sum to approach the analytical closed-form solution, 
Eq. (16). Convergence of the potential versus Zmax at 
these points are shown in Fig. [7j where the NFC are 
large, see below. Unlike previous approacheSj'^'^'^Sl the 
present method achieves a much better accuracy even at 
a lower /max. In contrast to Zhang' method, which hap- 
pen to produce fortuitously better potential for Imax = 4 
than Zmax = 6 near the corner of the cell (H-point), the 
overall quality of our potential improves consistently as 
/max is increased. Additionally, in all these other meth- 
ods, one needs to converge carefully the internal L'"*- 
sums; in most cases must be taken up to /^^x > •^^mlx^ 
and hence computationally expensive. However, Ham- 
merling et ahl^have shown that a multipole approach re- 
quires /"*x — S^max the van W. Morgan and Madelung 
models to achieve accuracy closer to our results. 

Again, the NFCs are the reason for a slower rate of 
convergence near the cell boundary, see Fig. |8j where the 
NFC contribution to the potential for an fee lattice are 
shown along the two symmetry directions with /max — 10. 
As before, this correction grows only beyond rmin £^nd 
get significant after tmt as the two densities in Eq. ([s]) 
are identical except outside the central cell where only 
p|f (r) 7^ 0. Unlike the van W. Morgan case, the NFC 



This extra distribution acts like an ionic atmosphere 
to screen the interactions between neighboring charges, 
which make these interactions now short-ranged, but all 
the Gaussian images must be summed to infinity. A 
closed-form solutiorl^l for the potential is given by 
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where Vq is an arbitrary constant and e is the Ewald 
parameter (controlling the width of the Gaussian in 
Eq. (14)), famously used to optimize the convergence of 
the sum used for screening, where part is done in real- 
space and part in k-space. Besides the on-site Gaussian, 
the erfc function requires summation over Gaussian tails 
contributing from neighboring sites, however many are 
non-zero. It can be verified that, with the constant of 
integration above, the potential is independent of e, as 
required, i.e., the first derivative with respect to e is zero. 

In Figure [6] we compare the numerical solution of the 
spatially-dependent potential from our general Eq. ([8| 
for /max = 0, 4, 6, 8, 10 to the numerical evaluation of the 



exact expression ( 16 ) for the jellium case for fee and bcc 
lattices. To assess the agreement, we used 15"^ Gauss 
points and 8 neighbor shells to evaluate Eq. (Isl). 
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FIG. 7. (Color online) rV(r) versus /max at high-symmetry 
points in eells of fee (top) and bee (bottom) for the Madelung 
jellium model. 
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FIG. 8. (Color online) For fee, the potential with (without) 
NFC along F— H and F— P for the Madelung jellium model. 
Other details are the same as in Fig. |4] 



along both the directions (especially along F-P) in the 
present case is relatively smaller, reflecting the distinct 
nature of the two models we have considered. 



Finally, we address the convergence properties of the 
Coulomb energy for the Madelung problem. By removing 
the self-energy arising in the blind application of Eq. ( 12 ) 



for the Madelung problem, a closed-form solution for the 
Coulomb energy U (for N unit cells) associated with the 
potential in Eq. (16) can be derived, i.e.. 



er/c(e|R„ 
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For convenience, rasa = (3rio/47r)^/^ is included, i.e., 
the radius for a sphere with equivalent unit cell volume 
flo, i.e., used in the atomic-sphere approximation (ASA). 
With this definition, U/{NZ'^r~}a) gives exactly 1.8 for 
the ASA Madelung problem, whereas the numerical eval- 
uation of Eq. ([T7| gives 1.79174723 (1.79185851) for fee 
(bcc), as found historically.^^ Using the potential and 
charge density within our Eqs. ([2|-([8|, we can evaluate 
the integrals for each VP and compare to the results of 
Eq. ([17]). 

Figure [9] shows the convergence of U versus Z^ax for 
fee and bcc lattices, compared to the exact result. For 
the Coulomb energy, the NFC do not have dramatic ef- 
fects, but there is error without them. No systematic 
error cancellation is possible, e.g., for Utcc-Ubcc, which is 
the well-known Ewald or "muffin-tin" corrections to the 
ASA structural energies. The Coulomb energy is almost 
correct by /max — 6 (error at 10"^ by Zmax = 8), and the 
convergence is monotonic, unlike when using multipole- 
based approaches with nested L sums, as shown by Ham- 
merling et al., where Z"*^ ^ 6/^*^ to achieve 10~^ ac- 
curacy comparable to our results without internal sums, 
which are very slowly convergent and numerically costly. 




max 



FIG. 9. (Color online) Coulomb energy for the Made lun g 
problem for fee and bcc, relative to the results from Eq. ( 17l. 



C. General Comments 

Our isoparametric integration avoids conditionally 
convergent summations, required in previous approaches, 
and provides a significantly more accurate and faster 
method for solving Poisson's equation, as detailed by 
the two cases. For molecular systems, a finite sum over 
atoms is required. For extended, solid-state systems, it 



also avoids FFTs, a limiting factor for large-atom cell 
calculations. In general, the present method is at least 
10(/'"' +_Lil-^VP times faster than any of the existing 
schemes .I^Em xhe factor (Z'"* -I- 1)^ comes from an ad- 
ditional internal L-sum (typically Z'"* ^ 6/°'^'), and the 
factor 10 is from use of isoparametric integration ver- 
sus shape functions, if used. In particular, for a system 
with Nvp sublattices, T''* ~ 8 - 10 will provide -lO-^Nyp 
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speed up. A direct comparison of CPU timings was de- 
tailed recentljK^ and shows that isoparametric integra- 
tion is 10^ faster and 10^ more accurate than that using 
shape functions. 

V. SUMMARY 

We have resolved the longstanding problem of an ac- 
curate, fast and efficient numerical solution of Poisson 
equation for electronic-structure codes with site-centered 
basis-sets. In particular, a proper calculation of the in- 
tercell potential has been developed that avoids trou- 
blesome multipole-type techniques that are conditionally 
convergent and we include accurately the correction term 
from the near cells, the so-called Near-Field Correction, 
where we have developed a physically intuitive and fast 
method to evaluate this correction also without multi- 
poles. The method provides machine-precision for po- 
tentials and Coulomb energy for systems described by 
arbitrarily-shaped, convex, space-filling VP, eliminates 
previous computational bottlenecks and convergence is- 
sues by employing isoparametric integration, scales as 
O(Nvp) and is easily parallelized. The method also 
avoids FFTs that do not scale well to very large cells. 



The method works for periodic solids, molecules (using 
extended VP) and materials containing imperfections or 
disorder. The general applicability and accuracy of the 
method was proved via two rigorous, analytic models 
that traverse from localized to extended densities. 



VI. ACKNOWLEDGEMENTS 

Research sponsored by the U.S. Department of Energy, 
Office of Basic Energy Science, Division of Materials Sci- 
ence and Engineering Division from contracts with DDJ 
(DEFG02-03ER46026) and seed funding with Ames Lab- 
oratory, which is operated for DOE by Iowa State Uni- 
versity under contract DE-AC02-07CH11358; from the 
"Center for Defect Physics" , an Energy Frontier Research 
Center, for DDJ to support a student who helped develop 
numerical integration method (Ref. llOp used here and in 
our EFRC's code. Work performed by BGW was under 
the auspices of the U.S. DOE by Lawrence Livermore Na- 
tional Laboratory under Contract DE-AC52-07NA27344. 
We also benefited from discussion with W.A. Shelton in 
our DOE/BES Computational Materials and Chemical 
Sciences Network, and D.M.C. Nicholson in the EFRC, 
to reproduce their method and results in Ref. |H1 



* emails: ddj,aftab@ameslab.gov,wilson9@llnl.gov| 

1 J. van W. Morgan, J. Phys. C: Solid State Kys. 10, 1181 
(1977). 

^ A. Gonis, Erik C. Sowa, and P. A. Sterne, Phys. Rev. 

Lett. 66, 2207 (1991). 
^ G. H. Schadler, Phys. Rev. B 45, 11314 (1992). 

* X.-G. Zhang, W. H. Butler, J. M. MacLaren, and J. van 
Ek, Phys. Rev. B 49, 13383 (1994). 

^ N. Stefanou, H. Akai and R. Zeller, Comput. Phys. Com- 
mun. 60, 231 (1990); Yang Wang, G.M. Stocks, and J.S. 
Faulkner,, Phys. Rev. B 49, 5028 (1994). 

^ M. Weinert, J. Math. Phys. 22, 2433 (1981); M. Weinert, 
et ai, Phys. Rev. 26, 4571 (1982). 

^ L. Vitos and J. KoUar, Phys. Rev. B 51, 4074 (1995). 

* D. M. C. Nicholson and W. A. Shelton, J. Phys.: Condens. 
Matter 14, 5601 (2002). 

® J. Zabloudil, R. Hammerling, L. Szunyogh, and P. Wein- 
berger, Electron Scattering m Solid Matter (Springer- 
Verlag, Berlin, 2005). 

° Aftab Alam, S. N. Khan, B. G. Wilson, and D. D. Johnson, 
Phys. Rev. B 84, 045105 (2011). 

^ John C. Slater, Insulators, Semiconductors and Metals, in 
Quantum Theory of Molecules and Solids, Vol. 3 (1967, 
McGraw-Hill, Inc., New York); see Chapters 4 and 9. 



^2 R. A. Sack, J. Math. Phys. 5, 260 (1964). 

" M. R. Pederson, D. V. Porezag, J. Kortus and D. C. Pat- 
ton, Phys. Status Solidi B 217, 197 (2000) [ NRLMO L: 
URL http : / /quantum. utep . edu/nrlmol/nrlm ol .html] . 
G. te Velde and E. J. Baere nds, Phy s. Rev. B 44, 7888 
(1991) [ADF: URL http: //www. scm. com/ . 
I. Dunlap, J. W. D. Connolly^ andT. R. Sabin, J. Chem. 
Phys. 71, 3396 (1979); ibid 71, 4993 (1979). 

^® Robert J. Buehler and Joseph O. Hirrschfelder, Phys. Rev. 
83, 3396 (1951); ibid 71, 149 (1951). 
J. Math Phys. 2, 825 (1961). 
Laurens Jansen, Phys. Rev. 110, 661 (1958). 

^® J. M. Dixon and R. Lacroix, J. Phys. A: Math, Nucl. Gen. 
6, 1119 (1973). 

2° W. I. van Rij, Phys. A: Math, Nucl. Gen. 8, 1164 (1973). 
R. A. Sack, SIAM J. Math. Anal. 5, 774 (1974). 
S. W. Leeuw, Proc. Roy Soc. A373, 27 (1980). 
P.P. Ewald, Ann. Phys. 64, 253 (1921). 
J.C. Slater and P. de Cicco, M.I.T. Quarterly Progress 
Report No. 50, Solid State and Molecular Theory Group, 
1963, p. 46. 

R. Hammerling, J. Zabloudil, L. Szunyogh. amd P. Wein- 
berger, Phil. Mag. 86(1), 25 (2006). 
'^^ e.g., Hans L. Skriver, Phys. Rev. B 31, 1909 (1985). 



